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Abstract 

We present a new Monte Carlo method which couples Path Integral for finite temperature protons with Quantum Monte Carlo 
for ground state electrons, and we apply it to metallic hydrogen for pressures beyond molecular dissociation. This method 
fills the gap between high temperature electron-proton Path Integral and ground state Diffusion Monte Carlo methods. Our 
data exhibit more structure and higher melting temperatures of the proton crystal than Car-Parrinello Molecular Dynamics 
results using LDA. We further discuss the quantum motion of the protons and the zero temperature limit. 
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1. Introduction 

Quantum Monte Carlo (QMC) methods have been 
developed for accurately solving the many-body 
Schrodinger equation. Zero temperature Variational 
Monte Carlo (VMC) and Diffusion Monte Carlo 
(DMC), and finite temperature Path Integral Monte 
Carlo (PIMC) are currently the most accurate and 
general methods for computing static properties of a 
quantum system [1,2]. They have been successfully 
applied to simple quantum many-body systems, such 
as the electron gas, hydrogen, and helium. 

Recently, there have been new attempts[3,4,5] 
to calculate properties of disordered systems such 
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as liquid hydrogen within QMC. For this system, 
VMC/DMC and PIMC are computationally too inef- 
ficient to provide definite answers, e.g. regarding the 
nature of the melting transition from liquid to solid or 
metal to insulator. Whereas PIMC calculations have 
been done at comparatively high temperatures [6], 
this method becomes computationally inefficient at 
temperatures lower than roughly 1/20 of the Fermi 
temperature. Zero temperature calculations (VMC, 
DMC) have been used for ground state calculations 
where both electronic and protonic degrees of free- 
dom are treated quantum mechanically [7,8]. How- 
ever, the convergence of these calculations suffers 
from the different masses of protons and electrons 
which introduce two time scales differing by three 
orders of magnitude, and, more important, low tem- 
perature properties are inaccessible by these ground 
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state methods. To fill this gap, the Coupled Electron- 
Ion Monte Carlo (CEIMC) has been developed [3,4] 
to combine a classical or quantum Monte Carlo sim- 
ulation of the nuclei at non-zero temperature with a 
QMC calculation for the electronic energies where the 
Born-Oppenheimer approximation helps to overcome 
the time scale problem. 

In Ref. [5], the CEIMC method has been applied 
to determine the equation of state of hydrogen for 
temperatures across the melting of the proton crystal. 
More structure and higher melting temperatures of the 
proton crystal compared to Car-Parrinello Molecular 
Dynamics (CPMD) results using LDA[9] have been 
found. In this paper, we shortly summarize the method 
and the results as reported in Ref. [5] and discuss in 
more detail the quantum effects of the protons [ 10, 1 1 ] . 



2. Method 

In the CEIMC method, the proton degrees of free- 
dom are advanced by a Metropolis algorithm in which 
the energy difference between the actual state S and 
the trial state 5" is computed by a Quantum Monte 
Carlo calculation The energies of the states are cal- 
culated within the Born-Oppenheimer approximation, 
where the electrons are assumed to remain in the 
ground state with respect to the actual protonic po- 
sitions. Since the Born-Oppenheimer energies E(S) 
and E(S') have to be sampled by a QMC calculation, 
they are affected by statistical noise which would bias 
the Monte Carlo sampling of the protons. At first sight 
one might expect that for an unbiased calculation one 
will need to reduce the accuracy of the energy differ- 
ence E(S) - E(S') much below k B T. However, it 
has been shown that unbiased sampling of the proton 
configurations can be efficiently achieved by using the 
penalty method[13], a generalization of the Metropo- 
lis algorithm, where detailed balance is satisfied on 
average. 

Since only differences of electronic energies are 
needed, we sample the electronic degrees of free- 
dom according to the sum of the electronic distribu- 
tion functions (e. g. the square of the trial wave func- 
tion in VMC) for the S and S' states, and we com- 
pute the energies for the two states as correlated sam- 
pling averages [3,4]. For the typical size of the proton 



moves (between O.OlAand 0.5Afor classical protons) 
this method is much more efficient than performing 
two independent electronic calculations [3,4]. 

An essential part of the CEIMC method is the choice 
of the trial wavefunction needed to calculate the Born- 
Oppenheimer energies. Variational Monte Carlo de- 
pends crucially on the quality of the trial wavefunc- 
tion. To go beyond VMC, we implemented a Reptation 
Quantum Monte Carlo algorithm (RQMC) [ 1 2] to sam- 
ple more accurately the electronic ground state. Simi- 
lar to DMC, RQMC projects the trial wavefunction on 
to the ground state within the Fixed-Node approxima- 
tion. A high quality trial wave functions is important 
to relax to the ground state with a very limited number 
of time slices and to provide accurate nodes. RQMC, 
being a Metropolis based method, is more easily used 
to compute energy differences; conversely, the corre- 
lated sampling method within DMC is more involved 
because of the branching step. 

To reduce finite size effects in metallic systems, 
we average over twisted boundary conditions (TABC) 
when computing electronic energies within CEIMC 
(i.e. we integrate over the Brillouin zone of the super 
cell) [15,4]. 

Quantum effects for protons are relevant at high 
pressure. We represent protons by imaginary time path 
integrals without considering the statistics of the pro- 
tons, (those effects are negligible in this temperature- 
density range.) For efficiency, it is important to mini- 
mize the number of protonic time slices. We have used 
the pair action of an effective proton-proton poten- 
tial and treated the difference between the true Born- 
Oppenheimer energy and the effective potential with 
the primitive approximation [2]. When coupled with 
TABC, rather than using all the k-points for each pro- 
tonic time slice, we can, randomly assign a subset of 
k-points to each protonic slice without introducing a 
detectable systematic effect. 



3. Results 



Comparison of CEIMC with CPMD We first con- 
sider classical protons. For classical protons it is pos- 
sible to compare the CEIMC results with previous 
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Fig. 1. Pair correlation function at r s = 1,T = 1000-ftT. 
Comparison between CEIMC-VMC-TABC with 32 protons, 
CEIMC-VMC-PBC with 54 protons and CPMD-LDA with 162 
protons (simulation with N p = 54 provides identical correlation). 
Data from CEIMC-VMC-TABC at T=2000K (stars) are also re- 
ported. 

Car-Parrinello molecular dynamics (CPMD) [9], the 
only difference is the method to calculate the poten- 
tial energy surface. Whereas in CEIMC the Born- 
Oppenheimer energies are calculated by QMC meth- 
ods, CPMD uses density functional theory (DFT) to 
calculate electronic energies. Both methods are in prin- 
ciple exact but rely on approximations of the unknown 
nodes of the trial wavefunctions in QMC and on the 
approximation of the unknown exchange-correlation 
energy functional in DFT. In Fig. 1 we compare the 
proton correlation function g pp (r) of both methods. 
The CEIMC results show more structure than CPMD 
of Ref. [9] using LDA. A better agreement is ob- 
served when CPMD results at temperature T are com- 
pared to CEIMC results at temperature 2T for 300 < 
T < 3000. Already previous studies suggested that 
the Born-Oppenheimer energy surface from DFT-LDA 
is flatter than the more accurate one from QMC. In 
Ref. [8] differences in energy among various crystal 
structures obtained within LDA were found smaller 
than DMC energies by roughly a factor of two. Fur- 
ther, in Ref. [3,4], the energies for random configu- 
rations of molecular hydrogen showed more variation 
within DMC than LDA. 

Note that the CEIMC curves of Fig. 1 are all 
obtained using VMC to obtain the electronic Born- 
Oppenheimer energies. RQMC requires roughly one 
order of magnitude more computer time then VMC. 
For this reasons large part of our results are based 
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Fig. 2. CEIMC-VMC-TABC with 54 protons. Protonic kinetic 
energy per particle at various densities versus temperature. The red 
line estimates the melting of the bcc crystal from the Lindemann 
ratios. 

on VMC calculation and RQMC is only exploited to 
estimate the systematic error of VMC. For the sys- 
tem with 54 protons at the V point, we have found 
no detectable differences in the correlation function 
between VMC and RQMC. 



The Quantum effects of the protons The quantum 
effects of the protons are summarized in Fig. 2 which 
shows the kinetic energy of the protons versus temper- 
ature for three different densities (r s — 1.2, 1.0,0.8) 
and the deviation from its classical value 3fceT/2. 
Furthermore, Fig. 3 compares g pp (r) for a classical 
and a quantum calculation which directly shows that 
the zero point motion is extremely important. The 
quantum effects not only affect the proton kinetic en- 
ergy but also increase the electronic energies and the 
configuration energy, at least at the variational level 
used in the present calculations. 



Trial wave functions and zero temperature results 
In this subsection we discuss the zero temperature 
limit by VMC/DMC calculations at T = of dy- 
namical protons. In most of our CEIMC calculations, 
analytic trial wave functions including backflow and 
three-body correlation [14] for the electronic degrees 
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Fig. 3. Comparision of gpp{r) between classical and quantum 
protons at r s = 1 and T = 500K (CEIMC-VMC-TABC). 

of freedom have been used. Here, we use this impor- 
tant ingredient of CEIMC to test their accuracy by 
comparing their energies for "dynamic" protons with 
the results of Natoli et al.[ll] in the ground state us- 
ing LDA trial wavefunctions. In these calculations, the 
LDA trial function was computed for a perfect bcc 
lattice and then modified for use within DMC calcu- 
lations of moving protons in order to avoid recalcula- 
tion of the LDA orbitals. In the present calculations, in 
addition to the backflow-threebody wavefunction for 
the electronic part of the wavefunction, the protonic 
part of the wavefunction contains a Jastrow correla- 
tion and has a Gaussian orbital localizing it to the bcc 
lattice sites (that was also in the Natoli calculation); 
exchange effects for protons due to statistics are ne- 
glected. In contrast to CEIMC, these DMC calcula- 
tions do not make use of the Born-Oppenheimer ap- 
proximation. From table 1 we see that the backflow 
wavefunctions have a lower energy by 2-3mH/atom 
within VMC and 2mH/atom within DMC. However, 
the analytical functions are particularly appropriate to 
our methods since their computational cost is much 
less then solving the Kohn-Sham equations for a ran- 
dom arrangement of 50 to 100 protons. 

To extend the data of the equation of state of Ref. [5] 
to the zero temperature limit, we performed VMC- 
TABC calculations at zero temperature. The results are 
summarized in table 2. The zero point energy is ob- 
tained by subtracting the energy in a static bcc lattice 
of protons from the energy of the dynamic electron- 
proton system. In the harmonic approximation, half of 
the zero point energy will contribute to the potential 



N 


wavefunction 


E v 


a 2 


Edmc 


16 


LDA 


-0.4678 (2) 








BF-A 


-0.4724(1) 


0.030 (2) 




54 


LDA 
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Table 1 



VMC and DMC Energies in h/atom (and variance per electron 
a) for dynamic protons and electrons at r s =1.31 at T = 0. LDA 
means LDA orbitals times an optimized 1 body factor and Jas- 
trow factor[8] used as trial wavefunction, BF-A are the analytical 
wavefunctions using backflow times an optimized gaussian for the 
protons. 
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0.18(1) 


1.0 
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20.02(6) 


0.023(1) 
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1.2 


-0.4637(4) 


0.873(2) 


-1.337(2) 


5.51(8) 


0.016(1) 


0.128(1) 



Table 2 



Energies for 54 dynamic protons and electrons at T = for 
various densities using VMC-TABC (averaged over 1000 twist 
angles); total (Etot), kinetic (E^in) and potential energies (E po t) 
are given in h/at, P is the pressure in Mbar, ZPE the zero point en- 
ergies and 7£ is the rms deviation divided by the nearest neighbor 
distance for a bcc lattice. 

and the other half to the kinetic energy of the pro- 
tons at zero temperature and we can compare with the 
protonic kinetic energies calculated in Ref. [5]. The 
zero temperature kinetic energies of the protons ob- 
tained is systematically higher than the zero temper- 
ature extrapolated results of Ref.[5]. It is known[7,8] 
that anharmonic effects are large for high pressure hy- 
drogen. The difference could also be due to limita- 
tions of the gaussian trial wavefunctions for the pro- 
tons used in the zero temperature VMC calculation. 
Indeed, in CEIMC, although electronic energies are 
calculated variationally, protons are represented by a 
path-integral; the CEIMC results are therefore more 
reliable. 
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